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We present an efficient implemention of a non-equilibrium Green function (NEGF) method for self- 
consistent calculations of electron transport and forces in nanostructured materials. The electronic 
structure is described at the level of density functional theory (DFT) using the projector augmented 
wave method (PAW) to describe the ionic cores and an atomic orbital basis set for the valence 
electrons. External bias and gate voltages are treated in a self-consistent manner and the Poisson 
equation with appropriate boundary conditions is solved in real space. Contour integration of the 
Green function and parallelization over k-points and real space makes the code highly efficient and 
applicable to systems containing several hundreds of atoms. The method is applied to a number 
of different systems demonstrating the effects of bias and gate voltages, multiterminal setups, non- 
equilibrium forces, and spin transport. 
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INTRODUCTION 

Electron transport across nanostructured interfaces 
is important in a range of different areas including 
nano-electronics, organic photovoltaics, and electroche- 
mistry. First-principles modelling of electron transport 
at the nano-scale has so far mostly been applied to mo- 
lecular junctions consisting of molecules contacted by 
metallic electrodes[l-5 . However, more recent applica- 
tions also include graphene nanoribbonsJSJ [7], semicon- 
ducting and metallic nanowires (5UTU]. and bulk tunnel- 
ling junctions for magneto-resistance and electrochemi- 
cal applications |TTJ [T2]. The rapid developments in these 
areas towards atomic-scale control of interface structures, 
and the continuing miniaturization of electronics com- 
ponents makes the development of efficient and flexible 
computational tools for the description of charge trans- 
port at the nano-scale an important endeavour. 

The vast majority of first-principles electron transport 
studies have been based on density functional theory 
(DFT) within the local density (LDA) or generalized gra- 
dient (GGA) approximations. This approach is in prin- 
ciple unjustified because the eigenvalues of the effecti- 
ve Kohn-Sham Hamiltonian do not represent the true 
quasiparticle energy levels. In particular, for tunneling 
junctions the energy gap between the highest occupied 
states and lowest unoccupied states is too small [T51 [T3] 
and this can lead to an overestimation of the conductan- 
ce. More accurate calculations incuding self-interaction 
corrections 1 15] and more recently the many-body GW 
approximation [TBHTH] yield conductance values in bet- 
ter agreement with experiments. On the other hand, the 
NEGF-DFT approach often provides a satisfactory quali- 
tative description [19] and its computational simplicity 
makes it a powerful tool for addressing non-equilibrium 
properties of complex systems. It should be mentioned 



that the formal problems associated with the use of DFT 
for transport are overcome by time-dependent DFT (TD- 
DFT) which allows for an, in principle, exact description 
of the (longitudinal) current due to an externally applied 
ficld|20j. However, it has been recently found that the 
standard TDDFT exchange-correlation potentials do not 
yield any improvement over the NEGF-DFT in terms of 
accuracy in predicting conductance|21|. 

In addition to the electronic current it is of interest 
to model the forces acting on the atoms under non- 
equilibrium conditions, i.e. under a finite bias voltage. 
Such forces ultimately determine the stability of current 
carrying molecular devices|22, 23J, but can also be exploi- 
ted to deliberately control the motion of single molecules 
by e.g. injecting electrons into the molcular orbitals using 
a scanning tunneling microscope (STM). 

In this paper we describe the implementation of the 
NEGF-DFT method in the GPAWpi] [55] electronic 
structure code. In GPAW the electronic states can be 
described either on a real space grid or using an ato- 
mic orbital basis set. For the NEGF calculations, the 
Green function is expanded in the atomic orbital basis 
while the Poisson equation is solved in real space. Con- 
tour integration and sparse matrix techniques together 
with parallelization over both k-points and real space is 
exploited for optimal efficiency. Although the basic ele- 
ments of our implementation are not new and have been 
described in earlier papers|26-28 , the possibility of ap- 
plying a general gate and/or finite bias voltage, the use 
of multiple leads, and inclusion of non-equilibrium forces 
on the ions provides a flexible and efficient computational 
platform for general purpose modelling of charge trans- 
port at the nano-scale and should be of interest to a large 
and growing community. 

This paper is organized as follows. In Section 2 the 
transport model and formalism are introduced. In Section 
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3 we describe the complex contour integration technique 
used to obtain the non-equilibrium electron density from 
the Green function. Section 4 describes the use of spar- 
se matrix methods, and in Section 5 we discuss the real 
space solution of the Poisson equation. A number of illu- 
strative applications are presented in Section 7. 
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Figur 1: A scattering region including the nanostructure of 
interest (e.g. a molecule) and part of the electrode atoms is 
sandwiched between two semi-infinite electrodes. Periodic bo- 
undary conditions are used in the x,y directions and open 
boundary conditions in the z direction. The electron poten- 
tial in the electrodes is periodic and can be obtained from 
a groundstate DFT calculation employing periodic boundary 
conditions in all directions. The Hartree potential at the scat- 
tering region boundary, which is used as boundary condition 
for the Poisson equation, is also obtained from the electrode 
calculation. The whole system can be subject to an external 
bias or gate voltage, and the electronic structure of the scat- 
tering region is calculated self-consistently in the presence of 
such external fields. 

The transport model is shown in Fig. [T] Following the 
standard approach, the system is divided into left and 
right electrodes and the central scattering region (see the 
detailed description in the caption). The Hamiltonian of 
the system is given by (all notation related to PAW met- 
hodology is consistent with earlier GPAW papers |2"4l |2"5] ) 
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where a denote the atoms in the system and z, j label the 
PAW projector functions of a given atom. Using a (non- 
orthogonal) atomic orbital basis set, the Hamiltonian can 
be written in the following generic form 
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The "on-site" Hamiltonian matrices of the electrodes, 
Hll (left) and Hrr (right), and the coupling matri- 
ces Hlc and Hue, can be obtained from a homoge- 
neous bulk calculation. If a bias voltage V is applied, 
the matrices corresponding to the left and right electro- 
des should be shifted by eV relative to each other, e.g. 
Hlc -> H L c + eVS L c and H L l -> H L L + eVS L L, where 
S denotes the overlap matrix. We assume that there is 
no coupling between basis functions belonging to diffe- 
rent electrodes. This assumption can be always satisfied 
by making the scattering region large enough. 
The retarded Green function is written as 



G r (s) = {(e)S - H CC - Z r L (e) - E£(e)] 



(3) 



The self-energies, ^ r L /R' represent the coupling to the 
electrodes and are obtained using the efficient decimation 
technique [29 . 

The lesser Green function is written as|30) 

G<(e) = G J -( £ )E < (e)G a ( £ ) + (1 + G r £ r )G^(l + E a G Q ). 

(4) 

The latter term is nonzero for truly bound states and 
vanishes for states acquiring any width. 

The pseudo density matrix (for the pseudo wave in the 
PAW framework) is the integral of G < 
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Here 9 = ksT and G is a contour for the integral to be 
discussed more in Sec. 3. 

The non-equilibrium electron density is obtained as 



n(r) 



(G) 



where <& v is an atomic orbital basis function and h a c is the 
atomic pseudo core density. As is standard in the PAW 
formalism a tilde indicates a smooth quantity as opposed 
to an all electron quantity. The smooth charge density is 
given by 



p(r)=n(r)+]T£Q? m ^ m (r), 



(7) 
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where Qf m are multipolc moments and <?" m (r) is a so- 
called shape function. The last term is the contribution 
to the charge density coming from the positively charged 
nuclei. 

The effective potential is found as 



5xc + E u °> 



(8) 



where the Coulomb potential satisfies the Poisson equa- 
tion V 2 u cou ; = — 47T/5, while v xc and v a are the exchange- 
correlation potential and zero-potential, respectively. v a 
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is a parameter chosen to smoothen v and which vanishes 
outside the augmentation sphere of atom q|31|. 

To obtain self-consistency we thus have the iteration 
process D — » p — » V e tt — > H — > D — > .... After conver- 
gence the current can be calculated by 
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where T L/R (e) = i(T, r L/R (e) ^ L 
of the current formula we refer the reader to Ref. 
(orthogonal basis) or Ref. [33J (non-orthogonal basis). 

The non-equilibrium force is obtained from the deriva- 
tive of the total energy with respect to atomic positions. 
In the PAW framework, the total energy is written 



E = E + ^{E a - E a ), 
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with 
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The force can be obtained as 
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We note that the expression given above does not include 
the recently discussed Berry phase contributions to the 
nonequilibrium force|23|. 



NUMERICAL DETAILS 



For the retarded Green function we use Gaussian 
quadrature by which a precision corresponding to a 
2N — 1 order polynomial can be obtained by N points. 
We use an adaptive method to find the energy points 
necessary to obtain a sufficient precision |3S]: for a given 
region [c— h, c + h], the integral Q of a function / can be 
estimated with the Gauss-Lobatto formula, 

Q = ^(f(c-h) + 5f(c- ~h) + 5f(c+^h) + f(c+h)) 

(13) 

and furthermore a Kronrod formula can be used to esti- 
mate the precision of the integral|36 
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Figur 2: The contour used for the Green function integral 
in the complex energy plane. The coordinates of the indi- 
cated points are: A(E m in, 0), where E m i n is less than all 
the eigen-energies of the system, which is usually taken as 
fiL — lOOeF since we only calculate the valence electrons sta- 
tes; T5(/-Ll + mksT, A), where m satisfies e _m w 0(a typical 
value of m is 10), ks is the Boltzmann constant, T is the 
electron temperature, A is normally between e n and e n +i, 
where et = (2i — l)n (with i a positive integer) is a pole 
of the Fermi-Dirac distribution function, thus the singulars 
below A should be counted when summing up the residues; 
C((iL—rnk B T, A); D(p R +mk B T, A); E(/i&— mfesT, r]), 77 is a 
infinitesimal to avoid the inversion divergence; F(fiR + mkBT, 

v). 



Contour integration technique 

The contour for the Green function integral in Eq. ^ 
is shown in Fig. [2] The retarded and lesser Green fun- 
ctions are integrated along the path AB (see Fig. [2]) in 
the upper half plane and the path EF closely above the 
real axis in the bias window|27] 134] . 



The difference between Q and Q' can be taken as the 
precision of the integral. 

The adaptive procedure to get the integral of the Green 
function in a region [c — h, c + h] is 

1. calculate Q and Q' , then compare the difference A 
with the tolerance S. 
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2. if A is smaller than 6, the integral is converged 
and Q is used as integral result. If not, divide 
the region [c — h,c + h] into three subregions 

[°- h ^ c - 7^]' [ c - 7f /l > c +7f /l ]> I c + 75 /l ' c + /l ] 
and redo the step 1 for each subregion until the 
integral is converged in the whole region. 

For the lesser Green function inside the bias window, 
we use the simple composite trapezoidal rule to obtain 
the integral. However, numerical errors can easily occur 
close to the real axis where the Green function has singu- 
larities. For this reason we apply the double contour met- 
hod introduced in Ref. 27: First, the integral of the re- 
tarded Green function is calculated along the path CD 
(Fig. [2| above the bias window, which is the spectrum 
for all the electronic states in the bias window S, then 
both electron spectrum G < and hole spectrum G > are 
integrated along the path EF, and we have S = D e + Dh 
according to the definition of the Green function, whe- 
re electron density D e and hole density Dh are obtained 
from the integral of G < and G > respectively. The numeric 
error, AS — S — D e — Dh is normally a non-zero quan- 
tity due to the integral insufficiency. As a correction, the 
error is distributed to D e and Dh by the matrix element 
weight. 

Sparse matrix handling 

Because the matrix inversion cost scales as TV 3 , where 
N is the dimension of the matrix, the matrix inversion 
turns out to be the main computational cost for large 
systems. Hence a sparse matrix method is implemented 
to obtain the Green function. 

We define a quenching layer as a slab whose left side 
has no overlap with the right side due to the finite cutoff 
in the range of the atomic orbitals. Hence an overlap or 
Hamiltonian matrix can be split into several blocks, with 
each block representing the onsite values of a quenching 
layer or the coupling between two adjacent quenching 
layers. Note that quenching layers here are different from 
principal layers used in the transport framework, where 
the latter is supposed to be repeatable as well. 

Physical quantities like density or transmission is often 
determined by fairly few blocks of a matrix. To see this 
consider the simple example of a two-probe system. In 
this case, the scattering region is divided into 5 quenching 
layers. Fig. [3] shows the sparse matrix structure of the 
overlap or the Hamiltonian matrix. The blocks outside 
the scattering region are from electrode calculations and 
always fixed. 

First we discuss how to obtain the real-space pseudo- 
density which can be obtained by a projection of the 
pseudo-density matrix as in Eq. [6] We see that if the sta- 
tes ^(r) and ^(r) have no overlap, the contribution 



from the pseudo-density matrix is zero, i.e., the white 
blocks in Fig. [3^i do not affect the density. So when we 
calculate the density matrix from the integral of the Gre- 
en' function, only the blue and green blocks in the Green 
function matrix (Fig. |3^) are necessary. There are real- 
ly two different parts, because two types of Green fun- 
ctions are involved when calculating the density matrix: 
the equilibrium part and the non-equilibrium part. We 
need the blue and green blocks of the retarded Green 
function for the former and of the Keldysh Green fun- 
ction for the latter. Through Eq. [4] and the finite extent 
of the self-energy matrix, which is only non-zero in the 
principle layers close to the electrodes, we see that the 
blue and green blocks in Fig. [3] in the Keldysh Green 
function matrix can be obtained from only the red blo- 
cks of the retarded Green function (Fig. [3j:). So when we 
do the matrix inversion to calculate the retarded Green 
function by Eq. |3j the red blocks in Fig. [3^i are neces- 
sary for energy points on the path EF in Fig. [2] and the 
blue and green blocks in Fig. [3b are necessary for energy 
points on the other path segments in Fig. [2] 

We can also see that the red and orange blocks in 
Fig. [3]3 are needed to calculate the density of states 
(DOS) by DOS(e) = -±Im(Tr(G(e+)S) and the pink 
blocks in Fig.|3]D are needed to calculate the transmission 
function T(e) = Tt\T L (e)G(E+)T R (e)G(e + )1]. 




Figur 3: Schematic of the matrix blocks. A, the Hamiltonian 
or overlap matrix, the blue and green blocks represent the 
on-site and coupling sub-matrices for the different quenching 
layers respectively; B, the Green function matrix when evalu- 
ating the density, DOS or transmission, the red and orange 
blocks represent the sub-matrices needed to calculate the den- 
sity matrix or DOS, and pink blocks are for the transmission 
coefficient; C, The red blocks in the retarded Green function 
matrix are necessary when calculating the Keldysh Green fun- 
ction. 

The formulas below provide a quick solution for the 
necessary blocks. Here we just consider this particular 
matrix(shown in Fig. [3| as an example to show how the 
method works. A general formalism, which works for ar- 
bitrary number of electrodes and arbitrary number of 
principal layers in each electrode, is introduced in Ref. 1371 

First, the central block N in Fig. |3]3 of the retarded 
Green function can be solved through the equations 

Q2 = ^22! Qi = (-^11 — L12Q2 L21) 1 (15) 

Q2 — Qi = {Ru ~ R12Q2 R21) 1 (16) 



5 



N = (M- JuQihi)- 1 , 

J=L,R 



(17) 



where and Rij are the blocks shown in Fig. [3^, rep- 
resenting the matrix eS — Hcc — ^l( £ ) — ^r( £ )- Then, 
for the remaining blocks of the retarded Green function 
matrix, we have to iterate the formulas 



G; . 
G 
G 



L 

i,3 
L 



— Qi (I — ^M-l^i-l.i) 

= (I ~ Gi.i-\L^_ li )Q^ 

= —QiLii^iGf_ 1 j 
= ,/., :,,()' <J < i) 



(18) 



where Gf is the block from the central part to electrode 



L and Gq is N in Eq. 
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the blocks from the central 
part to electrode R have~a similar solution. For all the 
required blocks, a quick solution can be obtained using 
a combination of the recursive formulas Eq. |18[ If we 
denote the number of quenching layers by n, the com- 
putational cost is roughly given by n times the cost of 
an inversion operation plus An times the cost of matrix 
multiplication|37 . 



Fixed boundary conditions 

The electronic potential of the metal electrodes will 
usually be very efficiently screened so that after only 
a few atomic layers into the electrodes we can assume 
that the potential is equal to the equilibrium potential 
plus/minus a possible constant bias potential. We shall 
apply open boundary conditions (in contrast to, say, pe- 
riodic ones) where the bias is applied by fixing the poten- 
tial values at the boundaries before solving the Poisson 
equation|38 . This procedure also allows for a net charge 
in the scattering region in which case the perturbation 
of the electron potential into the electrodes will of course 
be somewhat more long ranged. 

The Poisson equation V Vcoul — —Airp is solved in re- 
ciprocal space in the x and y directions while it is solved 
in real space in the ^-coordinate, i.e., in the the transport 
direction. Mathematically we have 



dz 2 



G 2 )v cou i(z,G) = -4irp{z,G), 



(19) 



where G is the vectors of the 2-d real grids used for the 
Fourier transformation. 

Eq. [19] is solved by the sparse matrix linear equation 
subroutine provided by the Lapack package. An advan- 
tage of this Poisson equation solution is the good paral- 
lelization behavior. The 2-d Fourier transformations are 
independent for the different real-space slices, and the 
linear equations Eq. [19] can be solved independently for 
different G-vectors. 
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Quantum capacitor 



3 



r t t 

■ H 

u ■ m * 

■ m 

M M M 



W T ■ 

■ KM 

■ • u ■ m 

M M M 
U M M 



i 




Figur 4: Upper panel: Capacitor model with the electrodes 
made of bcc sodium with the voltage drop along the (100) 
direction. The vacuum distance between the two electrodes is 
8 A, and the size of the unit cell in the transverse directions 
is 12. 7A x 12. 7A. The rectangle represents the scattering re- 
gion. Lower panel: The non-equilibrium part of the electron 
density (solid blue) and the induced effective potential (das- 
hed red) under a bias voltage of 5V, with the zero bias values 
as the reference. The calculated values are averaged over the 
transverse plane. 

We consider a simple capacitor system consisting of 
two semi-infinite Na electrodes separated by a vacuum 
gap, see upper part of Fig. [4] When a bias voltage is 
applied to the system, electrons accumulate/deplete on 
the Na(100) surfaces. According to the classic parallel 
plate capacitor model, the surface charge should be 



Q c i = eoVA/d, 



(20) 



where £q, A, d are the vacuum permittivity, area cross 
section of the unit cell, and the gap distance, respective- 
ly. We integrate the induced charge density in real spa- 
ce and obtain the net charge accumulation Q = 0.45e, 
which is close to the result of classic theory Q c i = 0.55e. 
One difference is that in the quantum theory, the charge 
decays from the surface into the bulk in an oscillatory 
fashion (Friedel oscillations), instead of being localized 
exactly on the surface as assumed in the classic model. 
We also note that at a distance of about 4 layers from the 
surface the values for both the potential and the charge 
are very close to their bulk values. Hence this calcula- 
tion confirms the screening approximation, namely that 
a few layers away from the surface or scattering region 
the potential has reached its bulk value. Finally, we note 
that the relatively high bias voltage of w hV is possible in 
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the present case where no current is flowing. On the other 
hand, the non-equilibrium states determining the current 
flow become increasingly difficult to calculate accurate- 
ly for larger bias values due to the insufficient integral 
of the Keldysh Green function in the bias window. As a 
consequence electron transport calculations are typically 
possible/reliable up to bias voltages of around 2 — 3V, 
depending on the transparency of the junction. 

Non-equilibrium forces 

The calculation of non-equilibrium forces is in prin- 
ciple a delicate problem involving non-conservative 
components [2U For highly conducting molecu- 

lar bridges an instability may occur which involves the 
Berry phase of the wave function. The description of such 
phenomena is beyond the scope of the usual NEGF+DFT 
framework, but in simpler cases, in particular in ca- 
ses with no or little current flow, the force expressions 
Eqs. |Tl|i"2"1 still apply. 

As an example, we show here a non-equilibrium for- 
ce calculation for a Au/N 2 /Au junction, where we can 
see the tendency towards molecular dissociation under a 
bias voltage. The structure (see upper part of Fig. [5]) is 
relaxed under zero bias until the maximum force is below 
O.OleV/A. When a positive bias is applied, electrons are 
redistributed over the molecule due to the electric field. 
Consequently, the two nitrogen atoms start to repel each 
other due to increased Coulomb repulsion which weakens 
the bond. The actual quantity of charge transfer to the 
molecule, which is about O.Ole for IV bias voltage on 
this system, shifts up the molecular energy spectrum, i.e. 
the energy levels follow the chemical potential of the left 
electrode (see Fig. [5] middle panel). The force is main- 
ly occurring only between the two nitrogen atoms, while 
there is no force induced between the electrode atoms 
and the N 2 molecule. Equivalently, a negative bias vol- 
tage shifts down the levels and pull electrons out of the 
N 2 molecule, and leads to an attractive force between the 
two nitrogen atoms. 

Electrostatic gate control 

One way of controlling the current flow through a nano- 
scale conductor is by electrostatic gating. This has be- 
en demonstrated experimentally for graphene, where a 
metal-insulator transition was induced by gating |4"0"ll41| . 
and for single-molecule junctions where the individual 
electronic levels were moved in energy relative to the 
Fermi level of the source/drain electrodes |421 14*5] . At the 
single molecule scale, the gate-molecule coupling is to a 
large extent determined by the device geometry with the 
screening effect playing an important role|44|. For nume- 
rical simulations, the typical method of applying a gate 
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Figur 5: The bias voltage effect on a nitrogen molecule be- 
tween two gold electrodes. Upper panel: the atomic structu- 
re of the system. The arrows represent the directions of the 
atomic forces generated by the bias voltage. Middle panel: 
PDOS(partial density of states) of the nitrogen molecule at 
0V and IV bias voltage. The Fermi level is located at zero 
and the red dashed lines show the location of the bias win- 
dow in the IV case. Lower panel: the magnitude of the non- 
equilibrium atomic forces as a function of bias voltage. 

is to add an external potential v g (r) to the effective po- 
tential of the system 

v(r)=v°(r) + v g (r). (21) 

We now consider the prototypical Au-BDT-Au jun- 
ction (see upper panel of Fig. [6| subject to three dif- 
ferent gate potentials, v g (z) (see middle panel of Fig.|6|. 
We note that the structure of the Au/BDT junction is 
presently being debated [45- 47 . However, for our purpose 
the simple model makes the sense. 

In the lower part of Fig. [6] we plot the resulting ef- 
fective gate potential Av(z) = v sc (z) — v® c {z) where the 
subscript sc denotes self-consistency, and the superscript 
denotes zero gate voltage. Due to the screening in the 
metal, the effective gate potentials only affect the mo- 
lecule region, and the narrower gate potential is seen to 
be less influenced by the self-consistency because it does 
not induce considerable charge transfer at the metal sur- 
faces - a charge transfer that otherwise tends to reduce 
the gate effect on the molecule. We note that the gate 
efficiency factor, a = Av(z)/v g (z), for these three poten- 
tials are about 0.8, 0.6, and 0.4 in the molecular region 
with the larger efficiency obtained for the more localized 
gate potential. The value of 0.4 obtained for the most 
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Figur 6: Effect of electrostatic gating of a benzene-1,4- 
dithiolate molecule between two gold fcc(lll) electrodes. Up- 
per panel: the atomic structure of the system. Middle panel: 
the applied external potential. Lower panel: the effective ad- 
ditional potential after self-consistency. 

delocalized potential is fairly close to an experimental 
studyf43| of the Au-BDT-Au system where an efficiency 
factor of 0.25 was reported. 

In the following we illustrate how the gate voltage can 
be used to tune the conductance of a molecular junction. 
It has recently been shown that the transport through 
the molecule anthraquinone is strongly suppressed due 
to destructive quantum interference occurring close to 
the Fermi level when the molecule is connected to me- 
tallic electrodes [48, [49]. The quantum interference leads 
to a dip in the transmission function inside the energy 
gap between the highest occupied (HOMO) and lowest 
unoccupied (LUMO) molecular orbitals. Hence a large 
on-off ratio is expected when shifting the molecular ener- 
gy levels by an external gate voltage. The upper panel of 
Fig. [7] shows the molecule connected to two gold fcc(lll) 
surfaces through Au-S bonds. The effective potential with 
the gate voltage -2V is shown in the middle panel. We see 
that the potential of the central part of the anthraquino- 
ne molecule is shifted less than the potential for the outer 
parts of the molecule. This is due to the fact that different 
parts of the molecule polarize differently as a consequen- 
ce of the detailed electronic structure. The HOMO is for 
example mainly localized at the connecting wires. The 
lower panel of Fig. [7] shows the change of transmission 
coefficient when a gate voltage of -2V is applied. Due to 
the characteristic interference dip in the transmission a 
large on-off ratio of about a factor of 1000 is achieved. 
Note that the relatively poor gate efficiency of around 0.1 
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Energy(eV) 



Figur 7: Gate-tuning the conductance of a molecular transi- 
stor. Upper panel: the Au-anthraquinone-Au structure. Mid- 
dle panel: the gate effect on the potential. Lower panel: the 
transmission coefficient at 0V and -2V gate voltages. 

is due to Fermi level pining of the HOMO. 

Spin transport 

In this section, we investigate the nonequilibrium- 
driven magnetic transition in the spin transport in 
zigzag graphene nanoribbon(ZGNR) which is proposed in 
Ref . |50| based on tight-binding calculations. The ZGNR's 
edge is spin-polarized and it has an anti-ferromagnetic 
spin configuration if its number of atomic layers is 
even[5T]. A gap about leV is opened between the diffe- 
rent spin states and makes the ZGNR a semiconductor. 
It was noticed by Denis et al. that the ZGNR's magnetic 
ordering is killed when the external bias voltage exceeds 
the size of the gap[50]. Here we reproduce this result for 
the graphene/ZGNR/graphene system shown in the up- 
per panel of Fig. [8] where a ZGNR(nn=8) is sandwiched 
between two semi-infinite graphene flakes. Under zero bi- 
as the PDOS of the central C atom along the ZGNR edge 
shows two peaks above and below the Fermi level, cor- 
responding to the different spin directions (middle panel 
in Fig. [8|. The distance between the two peaks is about 
1.6eV, and equals the band gap. When this bias voltage 
reaches l.OeV, the current starts to increase, see Fig. [9] 
At bias voltages 2.0eV the edge magnetic moment dis- 
appears very abruptly, and the current starts to increase 
even faster. 

Interestingly, in the tight-binding calculations presen- 
ted in Ref. |50j . both the magnetic moment and the 
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Figur 8: Spin transport in a zigzag graphene nanorib- 
bon(ZGNR) bridge connecting to two semi-infinite graphe- 
ne flakes. Upper panel: the atomic structure of the graphe- 
ne/ZGNR/graphene system. Middle panel: the PDOS of a C 
atom at the ribbon's zigzag edge under zero bias, the red and 
blue solid lines represent the spin-up and spin-down PDOS 
of the C atom at the center of the zigzag edge(marked with 
X), the red and blue dashed lines represent the spin-up and 
spin-down PDOS of the C atom next to the previous one 
at the zigzag edge(marked with Y), the green dashed line is 
the fermi level. Lower panel: the PDOS of a C atom at the 
ribbon's zigzag edge under bias voltage V = 2.4V. 



current show a very abrupt feature at the bias thres- 
hold, while in our calculation, the current increases rather 
smoothly. This can be explained by the non-equilibrium 
potential in the ab-initio calculation leads to a rehybri- 
dization and broadening of the spectral peaks, see lower 
panel of Fig. 8. We also note that in our calculation the 
disappearance of the magnetic moment is due to the two 
Stoner peaks moving into the bias window being half- 
occupied, different from the complete band collapse in the 
Ref. |Hn|)this is because our ZGNR is not long enough. 
We can see from the lower panel of Fig. [8] the gap shrinks 
more for the C atom further from the contact. 



Multi-terminal transport 

The expression for the Green function of the scattering 
region Eq. (§ can be straightforwardly extended to a 




Figur 9: Calculated current (blue squares) and edge magnetic 
moment per C atom (red circles) as a function of bias voltage 
for the ZGNR structure shown in Fig. [8] 



multi-terminal situation 



(22) 



where J is the index of the terminals. In contrast to the 
situation for a two-probe calculation, a zero boundary 
condition is applied for the effective potential for a multi- 
terminal system, and buffer atoms are used to represent 
semi-infinite leads. This approach to multiterminal trans- 
port has been previously investigated in Ref. [52) It should 
be noted that the self-energy of a lead has to be "rota- 
ted" by an orthogonal transformation when the lead is 
not along either the x, y or z axes. 

As an example, we consider a Cgo molecule connected 
to six linear carbon atomic chains. Fig.[lO|left) shows the 
projected DOS in real space evaluated at the Fermi level. 
The coverage suggests the scattering states are itinerant 
in the whole system and the contact between the carbon 
chain and the Cgo moelcule is transparent. Fig. 10 'right) 



shows a 2D averaged potential in a plane cutting through 
the Ceo molecule. 

A matrix indicating the transmission at the Fermi le- 
vel between the different leads is shown in the Fig. [TT] 
The matrix index notation represents the lead number- 
as shown in Fig. [Tp^left) . In particular, the diagonal cor- 
responds to back scattering, i.e. it gives the reflection 
probability. We can see that electrons are more easily 
transmitted between leads opposing each other, whereas 
the transmission decreases if the electron has to turn an 
angle during the scattering process. This intuitive pheno- 
menon can be explained by the quantum interference of 
the different partial waves. For the straightforward scat- 
tering, the quantum phases are the same for all the paths 
passing through the Ceo molecule, and the electron there- 
fore attains the greatest transmission probability[49j. 
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Figur 10: Left: The C6o-6-terminal structure and real space 
DOS at the Fermi level. Right: The averaged effective poten- 
tial projected onto a plane cutting through the C60 molecule. 



examples demonstrating electron transport under finite 
bias voltage, effect of electrostatic gating, spin transport, 
noncquilibrium forces, and multi-terminal leads. 
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Figur 11: The transmission function evaluated at the Fermi 
energy between the six leads. The element refers to the 

transmission from lead i to lead j; the diagonal element is the 
reflection coefficient. 



COMPUTATIONAL DETAILS 

For completeness we list the key input parameters and 
CPU timings for the examples presented in this paper in 
the table below. 



SUMMARY 

We have described the implementation of the NF- 
GF+DFT transport method in the GPAW code and illu- 
strated its performance through application to a number 
of different molecular junctions. The electronic structure 
is described within the PAW methodology which provi- 
des all-electron accuracy at a computational cost corre- 
sponding to pseudopotential calculations. The Green fun- 
ctions are represented in a basis set consisting of atomic- 
like orbitals while the Poisson equation with appropriate 
boundary conditions is solved in real space. The code 
is parallelized over k-points and real space domains and 
sparse matrix techniques are applied for maximal efficien- 
cy The flexibility of the method was illustrated through 
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System 


K-sampling 


Functional 


Basis 


Lead K-sampling 


processor number 


walltime(hour) 


Capacitor 


(4,4) 


LDA 


SZ 


(4,4,15) 


4 


0.5 


Au/N2/Au 


(2,2) 


PBE 


DZP(SZP) 


(2,2,15) 


32 


3 


Gate(BDT) 


(2,2) 


PBE 


DZP(SZP) 


(2,2,15) 


32 


6 


Gate(anthraquinon) 


(2,2) 


PBE 


DZP(SZP) 


(2,2,15) 


32 


8 


Graphene/ZGNR/Graphene 


(2,1) 


PBE 


SZP 


(2,2,7) 


32 


6 


C60 


(1,1) 


LDA 


SZ 


(1,1,15) 


32 


2 



Tabel I: Key parameters and CPU timings for the examples considered in this article. 
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